# --- Paths --------------------------------------------------------------------
base_dir <- Sys.getenv("PROJECT_DIR", unset = here::here())
rs_dir   <- file.path(base_dir, "Data", "Remote_Sensed")
shp_dir  <- file.path(base_dir, "Data", "Shape_Files")
out_dir  <- file.path(base_dir, "Processed_Data_h3_matrix_city")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

# --- Libraries ----------------------------------------------------------------
library(reshape2)
library(terra)
library(sf)
library(exactextractr)
options(stringsAsFactors = FALSE)

# --- Load rasters -------------------------------------------------------------
GHSb_2020 <- rast(file.path(rs_dir, "GHS_Built_100m_2020.tif"))
GHSp_2020 <- rast(file.path(rs_dir, "GHS_pop_100_2020.tif"))

# --- Load subplaces: CLEAN GPKG (same geometry as Section 4) ------------------
# This replaces the old raw SP_SA_2011.shp
sp <- terra::vect(file.path(shp_dir, "SP_SA_2011_clean.gpkg"))

# --- Load H3 grid shapefile (unchanged source) --------------------------------
h3 <- terra::vect(file.path(shp_dir, "UberH3_7.shp"))

# --- Assign CRS (assumed lon/lat) and project to raster CRS -------------------
# IMPORTANT: always overwrite CRS, do NOT try to be clever/conditional here.
# This is what the working version did and is robust to junk CRS metadata.
crs(sp) <- "EPSG:4326"
crs(h3) <- "EPSG:4326"

sp <- project(sp, GHSb_2020)
h3 <- project(h3, GHSb_2020)

# --- SP-level GHSL means + area ----------------------------------------------
sp_sf <- sf::st_as_sf(sp)

building_2020   <- exact_extract(GHSb_2020, sp_sf, "mean")
population_2020 <- exact_extract(GHSp_2020, sp_sf, "mean")

sp_df   <- as.data.frame(sp)
area_sp <- terra::expanse(sp)  # area in raster CRS units (as before)

sp_df <- cbind(
  sp_df,
  area_sp        = area_sp,
  building_2020  = building_2020,
  population_2020 = population_2020
)

# Keep only main variables and save
sp_keep <- c(
  "SP_CODE", "MP_CODE", "DC_MN_C", "SP_NAME", "MP_NAME",
  "DC_NAME", "PR_NAME", "building_2020", "population_2020", "area_sp"
)

sp_df1 <- sp_df[intersect(sp_keep, names(sp_df))]

saveRDS(sp_df1, file.path(out_dir, "sp_data.rds"))

# --- Per-city H3 × SP distance matrices --------------------------------------
city_vector <- c(
  "City of Cape Town",
  "City of Johannesburg",
  "eThekwini",
  "Nelson Mandela Bay",
  "City of Tshwane"
)

for (city in city_vector) {
  sp_city   <- subset(sp, sp$DC_NAME == city)
  sp_points <- centroids(sp_city)
  
  h3_city <- crop(h3, sp_city)
  h3_city <- terra::mask(h3_city, sp_city)
  
  # Force distance output to stay 2D even if one side = 1
  distance_h3_matrix_city <- matrix(
    distance(centroids(h3_city), sp_points),
    nrow = nrow(h3_city),
    ncol = nrow(sp_points)
  )
  
  rownames(distance_h3_matrix_city) <- as.vector(h3_city$hex7)
  colnames(distance_h3_matrix_city) <- as.vector(sp_city$SP_CODE)
  
  m <- melt(distance_h3_matrix_city)
  names(m) <- c("hex7", "SP_CODE", "Distance")
  
  out_base <- paste0("h3_sp_", gsub("\\s+", "_", city))
  saveRDS(m,  file.path(out_dir, paste0(out_base, ".rds")))
  write.csv(m, file.path(out_dir, paste0(out_base, ".csv")), row.names = FALSE)
}


